#include "RayTracer.h"
#include "tribox.h"
#include "material.h"
#include "numeric.h"
#include "optics.h"
#include <maya/MItDag.h>
#include <maya/MDagPath.h>

#include <boost/random/linear_congruential.hpp>
#include <boost/random/uniform_real.hpp>
#include <boost/random.hpp>
#include <boost/math/distributions.hpp>

#include <maya/MStatus.h>
#include <maya/MFnDagNode.h>
#include <maya/MFnCamera.h>
#include <maya/MFnLight.h>
#include <maya/MFnNurbsSurface.h>
#include <maya/MIOStream.h>
#include <maya/MPoint.h>
#include <maya/MVector.h>
#include <maya/MString.h>
#include <maya/MFnTransform.h>
#include <maya/MTransformationMatrix.h>
#include <maya/MMatrix.h>

#include <maya/MItDependencyNodes.h>
#include <maya/MLibrary.h>
#include <maya/MFileIO.h>

#include <maya/MFnMesh.h>
#include <maya/MFnMeshData.h>
#include <maya/MPointArray.h>
#include <maya/MItMeshVertex.h>
#include <maya/MItMeshPolygon.h>
//#include <maya/MFnPlugin.h>
#include <maya/MImage.h>
#include <maya/MFloatMatrix.h>
#include <maya/MCommonSystemUtils.h>

#include<maya/MFnAmbientLight.h>
#include<maya/MFnDirectionalLight.h>
#include<maya/MFnPointLight.h>
#include <maya/MPlug.h>
#include <maya/MPlugArray.h>
#include <maya/MMaterial.h>
#include <maya/MFnLambertShader.h>
#include <maya/MFnPhongShader.h>
#include <maya/MFnBlinnShader.h>
#include <maya/MTimer.h>
#include "Timer.h"

#include <vector>

const BBox BBoxSentinel = {std::numeric_limits<float>::max(), std::numeric_limits<float>::lowest(),
  std::numeric_limits<float>::max(), std::numeric_limits<float>::lowest(),
  std::numeric_limits<float>::max(), std::numeric_limits<float>::lowest()};

Visualizer RayTracer::s_viz;
double fRand(double fMin, double fMax)
{
    double f = (double)rand() / RAND_MAX;
    return fMin + f * (fMax - fMin);
}

float gen_random_float(float min, float max)
{
    boost::mt19937 rng;
    boost::uniform_real<float> u(min, max);
    boost::variate_generator<boost::mt19937&, boost::uniform_real<float> > gen(rng, u);
    return gen();
}

double inverseChi(int n, double beta)
{
	boost::math::chi_squared distrib(n);
	return boost::math::quantile(distrib, beta);
}

MString SST2STR(SuperSampleType sst){
  if(SuperSampleType_Unknown == sst)
    return MString("Unknown");
  else if (SuperSampleType_None == sst)
    return MString("none");
  else if (SuperSampleType_Uniform == sst)
    return SST_Uniform;
  else if (SuperSampleType_Jittered == sst)
    return SST_Jittered;
  else if (SuperSampleType_Random == sst)
    return SST_Random;
  else if (SuperSampleType_Adaptive == sst)
    return SST_Adaptive;
  else
    return MString("Unknown");
}


void* RayTracer::creator() 
{ 
  return new RayTracer; 
}

// thread job

// breaks down and creates the tasks

void DecomposeRayTrace(void *data, MThreadRootTask *root)
{
  TaskData* tdata = (TaskData*)data;  
  std::vector<ThreadData> tdvec;
  
  for (int i = 0; i < tdata->ranges.size(); i++)
  {
    ThreadData thr_data;
    thr_data.tidx = i;
    thr_data.chunk_range = tdata->ranges[i];
    thr_data.ray_tracer = tdata->ray_tracer;
    thr_data.pstat_data = &tdata->stats[i];
    tdvec.push_back(thr_data);
  }

  for(int i=0; i<tdata->thread_count; ++i)
  {
    MThreadPool::createTask(RayTraceChunk, (void*)&tdvec[i], root);
  }
  MThreadPool::executeAndJoin(root);
}

MThreadRetVal RayTraceChunk(void *data)
{
  ThreadData *thread_data = (ThreadData*)data;

  Range curr_range            = thread_data->chunk_range;
  
  RayTracer* ray_tracer       = thread_data->ray_tracer;
  RenderPlaneInfo& rp_info    = ray_tracer->m_rp_info;
  int arg_superSample         = ray_tracer->m_arg_superSample;
  VoxelSpace& vxspc           = ray_tracer->m_vxspc;
  std::vector<MeshInfo>& mesh_info_vec  = ray_tracer->m_mesh_info_vec;
  std::vector<LightInfo>& light_info_list = ray_tracer->m_light_info_list;
  int thread_count            = thread_data->ray_tracer->m_arg_thread_count;
  ThreadStatData* pstat_data  = thread_data->pstat_data;
  // timers
  
  unsigned char* pixel_array  = ray_tracer->m_img.pixels();
  int depth                   = ray_tracer->m_img.depth();

  MTimer group_timer;
  MTimer pixel_timer;
  LARGE_INTEGER freq;
  LARGE_INTEGER pixel_start;
  LARGE_INTEGER pixel_end;
  ::QueryPerformanceFrequency(&freq);
  
  int first_line = curr_range.y_start;
  
  MColor tolerance(ray_tracer->m_arg_tolerance,ray_tracer->m_arg_tolerance,ray_tracer->m_arg_tolerance);
  for(int y=curr_range.y_start; y<=curr_range.y_end; y += thread_count) 
  {
    if((y/thread_count)%10 == 0)
    {
      group_timer.endTimer();
      cout << "[Thread " << thread_data->tidx << "]: Processed lines " << first_line << " to " << y <<
        " (of " << curr_range.y_end << ") in " << group_timer.elapsedTime() << " s." << endl;
      cout.flush();
      first_line = y;
      group_timer.clear();
      group_timer.beginTimer();
    }
    
    for(int x=0; x<rp_info.ImgResX; ++x)
    {
      if(ray_tracer->m_arg_draw && ray_tracer->m_arg_drawOnly) {
        if(!((ray_tracer->m_arg_draw_x==x) && (ray_tracer->m_arg_draw_y==y))){
          continue;
        }
      }

      pixel_timer.clear();
      pixel_timer.beginTimer();
      //::QueryPerformanceCounter(&pixel_start);

      MColor avg_color;
      int color_sample_count = 0;

      // per-pixel super sampling loops
      float px_off = 1.0f/(arg_superSample+1); // pixel offset
      float d_x = (1.0f/(arg_superSample))*rp_info.PixelWidth; // distance in x between two samples
      float d_y = (1.0f/(arg_superSample))*rp_info.PixelHeight;

      //relevant if adaptive supersampling
      MColor incremental_sum_of_squares;
      MColor incremental_variance(FLT_MAX, FLT_MAX, FLT_MAX);
      int samples_number = 1;
      MColor incremental_mean;

      if(ray_tracer->m_arg_superSampleType == SuperSampleType_Adaptive)
      {
		  int choose = ray_tracer->m_arg_maxSR;
		  px_off = 1.0f/(choose+1);
		  d_x = (1.0f/(choose+1))*rp_info.PixelWidth;
		  d_y = (1.0f/(choose+1))*rp_info.PixelHeight;

        double x_offset1 = fRand(0,1);
        double y_offset1 = fRand(0,1);

        double perturb_x1 = fRand(-rp_info.PixelWidth/2, rp_info.PixelWidth/2);
        double perturb_y1 = fRand(-rp_info.PixelHeight/2, rp_info.PixelHeight/2);
		  //int rand_x1 = rand() % (choose+1) ;
		  //int rand_y1 = rand() % (choose+1) ;
		//double rand_x1 = double(choose+1)/2;
		//double rand_y1 = double(choose+1)/2;
        MVector pixelPos3d1 = rp_info.LeftBottom + 
          (double(y) + (0.5*rp_info.PixelHeight+perturb_y1)*px_off /*+perturb_y1*/)*(rp_info.PixelHeight*rp_info.upDirection) + 
          (double(x) + (0.5*rp_info.PixelWidth+perturb_x1)*px_off /*+perturb_x1*/)*rp_info.PixelWidth*rp_info.rightDirection;

        /*MVector pixelPos3d1 = rp_info.LeftBottom + 
        (double(y) + y_offset1)*(rp_info.PixelHeight*rp_info.upDirection) + 
        (double(x) + x_offset1)*rp_info.PixelWidth*rp_info.rightDirection;*/
        MFloatPoint raySource1;
        raySource1.setCast(pixelPos3d1);
        MVector rayDir1;
        if(rp_info.isOrtho){
          rayDir1 = rp_info.viewDirection;
          rayDir1.normalize();}
        else{
          MVector sampleRay = pixelPos3d1 - rp_info.eyePoint;
          rayDir1 = MFloatVector(sampleRay);
          rayDir1.normalize();}

        // call raytrace
        int traceDepth1=0;
        MColor out_color1;
        bool draw1 = false;
        if(ray_tracer->m_arg_draw){
          if((ray_tracer->m_arg_draw_x==x) && (ray_tracer->m_arg_draw_y==y)){
            draw1 = true;
            cout << "[Thread " << thread_data->tidx << "]: Drawing for pixel (" << x << ", " << y << ")" << endl;
          }
        }

        CIStats ci_stats_local;
        int ray_count_local;
        std::vector<int> sample_ray_depths_local;
        sample_ray_depths_local.clear();
        bool rt_success = ray_tracer->RayTrace(raySource1, rayDir1, traceDepth1, 1.0f, out_color1, 
          ci_stats_local, ray_count_local, sample_ray_depths_local, draw1, NULL);
        pstat_data->ci_stats += ci_stats_local;
        pstat_data->total_rays += ray_count_local;
        for(int s=0; s<sample_ray_depths_local.size(); ++s){
          pstat_data->sample_ray_depths.push_back(sample_ray_depths_local[s]);
        }
        incremental_mean = out_color1;
		if(out_color1.b!=0)
			bool stop = true;

        while((incremental_variance.r >= tolerance.r*ray_tracer->m_chiInv_R[int((ray_tracer->beta_R-1.0/double(CHI_INV_COUNT))*CHI_INV_COUNT)][samples_number] 
        || incremental_variance.g >= tolerance.g*ray_tracer->m_chiInv_R[int((ray_tracer->beta_G-1.0/double(CHI_INV_COUNT))*CHI_INV_COUNT)][samples_number]
        || incremental_variance.b >= tolerance.b*ray_tracer->m_chiInv_R[int((ray_tracer->beta_B-1.0/double(CHI_INV_COUNT))*CHI_INV_COUNT)][samples_number] 
        || samples_number < ray_tracer->m_arg_minSR) && samples_number <= ray_tracer->m_arg_maxSR )
        { 
			if(samples_number==1){
				incremental_variance.r=0;
				incremental_variance.g=0;
				incremental_variance.b=0;
			}
          double x_offset = fRand(0,1);
          double y_offset = fRand(0,1);
          double perturb_x = fRand(-rp_info.PixelWidth/2, rp_info.PixelWidth/2);
        double perturb_y = fRand(-rp_info.PixelHeight/2, rp_info.PixelHeight/2);
		  int rand_x = rand() % (choose+1) ;
		  int rand_y = rand() % (choose+1) ;
		  //double rand_x = double(choose+1)/2;
		//double rand_y = double(choose+1)/2;
         MVector pixelPos3d = rp_info.LeftBottom + 
          (double(y) + (0.5*rp_info.PixelHeight+perturb_y)*px_off /*+perturb_y1*/)*(rp_info.PixelHeight*rp_info.upDirection) + 
          (double(x) + (0.5*rp_info.PixelWidth+perturb_x)*px_off /*+perturb_x1*/)*rp_info.PixelWidth*rp_info.rightDirection;

          /*MVector pixelPos3d = rp_info.LeftBottom + 
          (double(y) + y_offset)*(rp_info.PixelHeight*rp_info.upDirection) + 
          (double(x) + x_offset)*rp_info.PixelWidth*rp_info.rightDirection;*/

          MFloatPoint raySource;
          raySource.setCast(pixelPos3d);

          MVector rayDir;
          if(rp_info.isOrtho)
          {
            rayDir = rp_info.viewDirection;
            rayDir.normalize();
          }
          else
          {
            MVector sampleRay = pixelPos3d - rp_info.eyePoint;
            rayDir = MFloatVector(sampleRay);
            rayDir.normalize();
          }

          // call raytrace
          int traceDepth=0;
          MColor out_color;
          samples_number++;

          bool draw = false;
          if(ray_tracer->m_arg_draw){
            if((ray_tracer->m_arg_draw_x==x) && (ray_tracer->m_arg_draw_y==y)){
              //if((91==x) && (205==y)){
              draw = true;
              cout << "[Thread " << thread_data->tidx << "]: Drawing for pixel (" << x << ", " << y << ")" << endl;
            }
          }
          
          int found_poly_intersections, total_vox_rays, total_poly_rays;
          sample_ray_depths_local.clear();
          bool rt_success = ray_tracer->RayTrace(raySource, rayDir, traceDepth, 1.0f, out_color, 
            ci_stats_local, ray_count_local, sample_ray_depths_local, draw, NULL);
          pstat_data->ci_stats += ci_stats_local;
          pstat_data->total_rays += ray_count_local;
          for(int s=0; s<sample_ray_depths_local.size(); ++s){
            pstat_data->sample_ray_depths.push_back(sample_ray_depths_local[s]);
          }

          MColor D = out_color - incremental_mean;
          incremental_mean += (D/float(samples_number));
          incremental_sum_of_squares += D*(out_color - incremental_mean);
          incremental_variance = incremental_sum_of_squares/float(samples_number);
        }
		pstat_data->total_samples += samples_number;
		pstat_data->pixel_samples.push_back(samples_number);
      }
      else
      {
        for(int sx=0; sx<arg_superSample; ++sx)      
        {
          for(int sy=0; sy<arg_superSample; ++sy)
          {
            // set up the source point and direction
            MVector pixelPos3d;
            if(ray_tracer->m_arg_superSampleType == SuperSampleType_Uniform || ray_tracer->m_arg_superSampleType == SuperSampleType_None )
            {
              pixelPos3d = rp_info.LeftBottom + 
              (double(y) + (sy+1)*px_off)*(rp_info.PixelHeight*rp_info.upDirection) + 
              (double(x) + (sx+1)*px_off)*rp_info.PixelWidth*rp_info.rightDirection;
            }
            else
            {
              if(ray_tracer->m_arg_superSampleType == SuperSampleType_Jittered)
              {
                double perturb_x = fRand(-d_x/2, d_x/2);
                double perturb_y = fRand(-d_y/2, d_y/2);
                pixelPos3d = rp_info.LeftBottom + 
                  (double(y) + (double(sy)+1.0)*px_off +perturb_y)*(rp_info.PixelHeight*rp_info.upDirection) + 
                  (double(x) + (double(sx)+1.0)*px_off+perturb_x )*rp_info.PixelWidth*rp_info.rightDirection;
              }
              else
                if(ray_tracer->m_arg_superSampleType == SuperSampleType_Random)
                {
                  double x_offset = fRand(0,rp_info.PixelWidth);
                  double y_offset = fRand(0,rp_info.PixelHeight);
                  pixelPos3d = rp_info.LeftBottom + 
                    (double(y) + (x_offset+1.0))*(rp_info.PixelHeight*rp_info.upDirection) + 
                    (double(x) + (y_offset+1.0))*rp_info.PixelWidth*rp_info.rightDirection;
                }
            }

            MFloatPoint raySource;
            raySource.setCast(pixelPos3d);

            MVector rayDir;
            if(rp_info.isOrtho)
            {
              rayDir = rp_info.viewDirection;
              rayDir.normalize();
            }
            else
            {
              MVector sampleRay = pixelPos3d - rp_info.eyePoint;
              rayDir = MFloatVector(sampleRay);
              rayDir.normalize();
            }

            // call raytrace
            int traceDepth=0;
            MColor out_color;
            bool draw = false;
            if(ray_tracer->m_arg_draw){
              if((ray_tracer->m_arg_draw_x==x) && (ray_tracer->m_arg_draw_y==y)){
                //if((91==x) && (205==y)){
                draw = true;
                cout << "[Thread " << thread_data->tidx << "]: Drawing for pixel (" << x << ", " << y << ")" << endl;
              }
            }

            int found_poly_intersections, total_vox_rays, total_poly_rays;
            CIStats ci_stats_local;
            int ray_count_local;
            std::vector<int> sample_ray_depths_local;
            bool rt_success = ray_tracer->RayTrace(raySource, rayDir, traceDepth, 1.0f, out_color, 
              ci_stats_local, ray_count_local, sample_ray_depths_local, draw, NULL);
            pstat_data->ci_stats += ci_stats_local;
            pstat_data->total_rays += ray_count_local;
            avg_color += out_color;
            for(int s=0; s<sample_ray_depths_local.size(); ++s){
              pstat_data->sample_ray_depths.push_back(sample_ray_depths_local[s]);
            }

            ++color_sample_count;
          } // for sy
        } // for sx
		pstat_data->total_samples += (arg_superSample*arg_superSample);
		pstat_data->pixel_samples.push_back(arg_superSample*arg_superSample);
      }
      // we have the sum; now divide

      if(ray_tracer->m_arg_superSampleType == SuperSampleType_Adaptive)
      {
        pixel_array[y*rp_info.ImgResX*depth + x*depth + 0] = RoundToInt(incremental_mean.r*255);
        pixel_array[y*rp_info.ImgResX*depth + x*depth + 1] = RoundToInt(incremental_mean.g*255);
        pixel_array[y*rp_info.ImgResX*depth + x*depth + 2] = RoundToInt(incremental_mean.b*255);
        pixel_array[y*rp_info.ImgResX*depth + x*depth + 3] = 255;
      }
      else
      {
        avg_color /= float(color_sample_count);
        pixel_array[y*rp_info.ImgResX*depth + x*depth + 0] = RoundToInt(avg_color.r*255);
        pixel_array[y*rp_info.ImgResX*depth + x*depth + 1] = RoundToInt(avg_color.g*255);
        pixel_array[y*rp_info.ImgResX*depth + x*depth + 2] = RoundToInt(avg_color.b*255);
        pixel_array[y*rp_info.ImgResX*depth + x*depth + 3] = 255;
      }
      
      pixel_timer.endTimer();
      //::QueryPerformanceCounter(&pixel_end);
      //float elapsed_time = float(double(pixel_end.QuadPart-pixel_start.QuadPart) / double(freq.QuadPart));
      float elapsed_time = pixel_timer.elapsedTime();
      pstat_data->pixel_times.push_back(elapsed_time);

    } // for x->ImgResX

  } // for y->ImgResY



  return (MThreadRetVal)0;
}

RayTracer::RayTracer()
{
  m_chiInv_R = NULL;
}

//---------------------------------------------------------------------------
MStatus RayTracer::doIt(const MArgList& argList) 
{
  // create the output directory if it doesn't already exist
  MString img_dir = buildOutputDirectoryName();
  MCommonSystemUtils::makeDirectory(img_dir);
  int charlen = img_dir.length();
  const char* tmp = img_dir.asChar(charlen);

  SYSTEM_INFO sys_info;
  ::GetSystemInfo(&sys_info);
  m_core_count = sys_info.dwNumberOfProcessors;


	std::fstream file;

	std::string str(tmp);
	str += "performance_report.txt";
	file.open(str,ios::out);

  int ret_width;
  int ret_height;
  int ret_voxelFactor;
  int ret_superSample;
  SuperSampleType ret_superSampleType;
  int ret_samplingRate;
  float ret_tolerance;
  int ret_minSR;
  int ret_maxSR;
  int ret_rayDepth;
  int ret_verbose;
  int ret_draw_x;
  int ret_draw_y;
  bool ret_drawOnly;
  int ret_thread_count;
  
  parseArgs(argList, ret_width, ret_height, ret_voxelFactor, 
    ret_superSample, ret_superSampleType, ret_samplingRate, 
    ret_tolerance, ret_minSR, ret_maxSR, ret_rayDepth, ret_verbose, 
    ret_draw_x, ret_draw_y, ret_drawOnly, ret_thread_count);

  m_arg_width           = (ret_width != -1) ? ret_width : 1920;
  m_arg_height          = (ret_height != -1) ? ret_height : 1080;
  m_arg_voxelFactor     = (ret_voxelFactor != -1) ? ret_voxelFactor : 1;
  m_arg_superSample     = (ret_superSample != -1) ? ret_superSample : 1;
  m_arg_superSampleType = (ret_superSampleType != SuperSampleType_None) ? ret_superSampleType : SuperSampleType_Uniform;
  m_arg_samplingRate    = (ret_samplingRate != -1) ? ret_samplingRate : 1;
  m_arg_tolerance       = (ret_tolerance > 0) ? ret_tolerance : 1;//1e-6f; // default based on experiments
  m_arg_minSR           = (ret_minSR != -1) ? ret_minSR : 1;
  m_arg_maxSR           = (ret_maxSR != -1) ? ret_maxSR : 128;
  m_arg_rayDepth        = (ret_rayDepth != -1) ? ret_rayDepth : 5;
  m_arg_verbose         = (ret_verbose != -1) ? ret_verbose : 0;

  m_arg_draw_x              = (ret_draw_x != -1) ? ret_draw_x : -1;
  m_arg_draw_y              = (ret_draw_y != -1) ? ret_draw_y : -1;
  m_arg_draw                = (m_arg_draw_x != -1) && (m_arg_draw_y != -1);
  m_arg_drawOnly        = ret_drawOnly;
  m_arg_thread_count    = (ret_thread_count > 0) ? ret_thread_count : m_core_count;

  cout << endl << endl << endl << endl;
  cout << "==============================" << endl;
  cout << "raytrace" << endl;
  cout << "   width: " << m_arg_width << endl;
  cout << "   height: " << m_arg_height << endl;
  cout << "   voxel factor: " << m_arg_voxelFactor << endl;
  cout << "   super sampling: " << m_arg_superSample << endl;
  cout << "   super sample type: " << SST2STR(m_arg_superSampleType).asChar() << endl;
  cout << "   sampling rate: " << m_arg_samplingRate << endl;
  cout << "   tolerance: " << m_arg_tolerance << endl;
  cout << "   min/max SR: " << m_arg_minSR << "/" << m_arg_maxSR << endl;
  cout << "   ray depth: " << m_arg_rayDepth << endl;
  cout << "   verbose: " << (m_arg_verbose?"on":"off") << endl;
  cout << "   draw: " << (m_arg_draw?"yes":"no") << "(" << m_arg_draw_x << "," << m_arg_draw_y << ")" << endl;
  cout << "   threads: " << m_arg_thread_count << endl;
  SYSTEMTIME sys_time;
  ::GetLocalTime(&sys_time);
  cout << "Called on " << sys_time.wYear << "/" << sys_time.wMonth << "/" << sys_time.wDay <<  " at " <<
    sys_time.wHour << ":" << sys_time.wMinute << ":" << sys_time.wSecond << "." << sys_time.wMilliseconds << endl;

  computeChiInvForBetas();
  beta_R = 0.01; 
  beta_G = 0.01;
  beta_B = 0.01;
  averageSamplingRate = 0;

  // total timer
  MTimer total_timer;
  total_timer.clear();
  total_timer.beginTimer();

  // init timer
  double init_time = 0.0;
  Timer init_timer;
  init_timer.clear();
  init_timer.beginTimer();
  file << "Preprocessing begun.." << endl;
  file << "Triangulating the meshes in the scene.." << endl;

  // triangulate meshes
  triangulate_meshes();

  init_timer.endTimer();
  file << "Elapsed time: " << init_timer.elapsedTime() << endl;
  init_time += init_timer.elapsedTime();
  file << "Reading meshes information.." << endl;
  init_timer.clear();
  init_timer.beginTimer();

  // create a vector of structures with MeshInfo
  
  create_meshinfo_list(m_mesh_info_vec);

  // create a vector of structures with triangle info
	int total_number_of_triangles = 0;  
  create_triangles_list(m_triangles_vector, total_number_of_triangles, (bool)m_arg_verbose);

  // get a bounding box enclosing all meshes in the scene
  float scale_factor = 1.1f;
  BBox final_bbox;
  get_scene_bounding_box(final_bbox, scale_factor);

  init_timer.endTimer();
  init_time += init_timer.elapsedTime();
  file << "Done reading meshes information!" << endl;
  file << "Elapsed time: " << init_timer.elapsedTime() << endl;
  file << "Voxel space initialization begun.." << endl;
  init_timer.clear();
  init_timer.beginTimer();

  // create the voxel space  
  //VoxelSpace vxspc;
  m_vxspc.Init(final_bbox, m_arg_voxelFactor);
  if(!m_vxspc.AssignTrianglesToVoxels(m_triangles_vector))
  {
    MString msg("Could not Voxelize. VoxelSpace probably not initialized.");
    MGlobal::displayInfo(msg);
    file.close();

    return MS::kFailure;
  }

  init_timer.endTimer();
  init_time += init_timer.elapsedTime();
  file << "Voxel space initialization ended!" << endl;
  file << "Elapsed time: " << init_timer.elapsedTime() << endl;
  init_timer.clear();
  init_timer.endTimer();
  file << "Reading Camera information.." << endl;
  init_timer.beginTimer();

  // get render plane info  
  MStatus rpStat = get_renderplane_info(m_arg_width, m_arg_height, m_rp_info);
  if(MS::kFailure == rpStat)
  {
    // probably no cameras
    MString msg("Could not create render plane. Probably no non-default cameras.");
    MGlobal::displayInfo(msg);
    file.close();
    return MS::kFailure;
  }

  init_timer.endTimer();
  init_time += init_timer.elapsedTime();
  file << "Done reading camera information!" << endl;
  file << "Elapsed time: " << init_timer.elapsedTime() << endl;
  file << "Reading lights information.." << endl;
  init_timer.clear();
  init_timer.beginTimer();

  // get all lighting information. For each light:
  create_lightinfo_list(m_light_info_list);

  init_timer.endTimer();
  init_time += init_timer.elapsedTime();
  file << "Reading lights information ended!" << endl;
  file << "Elapsed time: " << init_timer.elapsedTime() << endl;
  file << "Ray Casting started" << endl;
  init_timer.clear();
  init_timer.beginTimer();

  // create the image
  MStatus stat_img;
  //MImage  img;
  stat_img                    = m_img.create(m_rp_info.ImgResX, m_rp_info.ImgResY);
  int depth                   = m_img.depth();
  unsigned char* pixel_array  = m_img.pixels();
  memset(pixel_array, 0, m_rp_info.ImgResX*m_rp_info.ImgResY*depth);

  init_timer.endTimer();
  init_time += init_timer.elapsedTime();
  file << "Image creation ended!" << endl;
  file << "Elapsed time: " << init_timer.elapsedTime() << endl;


  file << "Threading setup started" << endl;
  init_timer.clear();
  init_timer.beginTimer();

  int thread_count = m_arg_thread_count;
  float frows_per_thread = float(m_rp_info.ImgResY)/thread_count;
  int rows_per_thread   = RoundToInt(ceil(frows_per_thread));
  std::vector<Range> row_ranges;
  for(int i=0; i< thread_count; ++i)
  {
    int ystart = i;
    int yend = m_rp_info.ImgResY-1;
    Range new_range(ystart, yend);
    row_ranges.push_back(new_range);
  }

  // thread-ify
  MStatus t_stat = MThreadPool::init();
  if( MStatus::kSuccess != t_stat ) {
    MString str = MString("Error creating threadpool");
    MGlobal::displayError(str);
    file.close();
    return MStatus::kFailure;
  }

  // build the data for the task
  TaskData task_data(this, row_ranges);

  init_timer.endTimer();
  init_time += init_timer.elapsedTime();
  file << "Threading setup ended!" << endl;
  file << "Elapsed time: " << init_timer.elapsedTime() << endl;
  file << "Preprocessing ended!" << endl;
  file << "Total time for preprocessing " << init_time << endl;
  init_timer.clear();

  file << "Ray Tracing started" << endl;

  MTimer raytrace_timer;
	raytrace_timer.beginTimer(); 

  // create the parallel region
  MThreadPool::newParallelRegion(DecomposeRayTrace, (void*)&task_data);

  raytrace_timer.endTimer();
  double raytrace_elapsed = raytrace_timer.elapsedTime();
  file << "RayTracing ended!!" << endl;

  MThreadPool::release();
  MThreadPool::release();

  // finally, write the image
  MString output_filename = img_dir + MString("scene.tif");
  stat_img = m_img.writeToFile(output_filename,"tif");
  if(stat_img) {
    MString msg = MString("Successfully wrote scene file ")+output_filename;
    MGlobal::displayInfo(msg);
    cout << msg.asChar() << endl;
  } else {
    MGlobal::displayInfo(MString("Error writing scene file ")+output_filename);
  }

  // release the image
  m_img.release();
  for(int m=0; m<m_mesh_info_vec.size(); ++m)
  {
    if(m_mesh_info_vec[m].m_material.m_isTextured){
      m_mesh_info_vec[m].m_material.release_image();
    }
  }


  // combine thread stat data
  std::vector<float> all_pixel_times;
  std::vector<int> all_pixel_samples;
  std::vector<int> all_path_lengths;
  int total_samples = 0;
  int total_rays = 0;
  int total_depths = 0;
  
  CIStats cis_stat_all;  for(int i=0; i<task_data.thread_count; ++i){
	  total_samples += task_data.stats[i].total_samples;
    for(int t=0; t<task_data.stats[i].pixel_times.size(); ++t){
      all_pixel_times.push_back(task_data.stats[i].pixel_times[t]);
	    all_pixel_samples.push_back(task_data.stats[i].pixel_samples[t]);
      all_path_lengths.push_back(task_data.stats[i].sample_ray_depths[t]);
      total_depths += task_data.stats[i].sample_ray_depths[t];
    }
    total_rays += task_data.stats[i].total_rays;
    cis_stat_all += task_data.stats[i].ci_stats;
  }
  cout << "Count of all_pixel_times: " << all_pixel_times.size() << endl;
  cout << "Total rays: " << total_rays << endl;
  cout << "Total ray-polygon intersection tests: " << cis_stat_all.total_ray_poly_tests << endl;
  cout << "Total voxels traversed: " << cis_stat_all.total_voxels_traversed << endl;
  cout << "Total positive intersection tests: " << cis_stat_all.found_poly_intersections << endl;
  cout << "Total samples: " << total_samples << "." << endl;
  cout << "Average Ray Depth: " << double(total_depths) / double(total_samples) << endl;

  int total_pixels = (m_rp_info.ImgResX*m_rp_info.ImgResY);
  double time_per_pixel_overall = raytrace_elapsed / total_pixels;
  cout << "Total pixels in image: " << total_pixels << endl;
  cout << "Time per pixel (end_time-start_time / total_pixels): " << time_per_pixel_overall << endl;

  float total_pixel_time = 0.0f;
  for(int i=0; i<all_pixel_times.size(); ++i){
    total_pixel_time += all_pixel_times[i];
  }
  float avg_pixel_time_approx = (total_pixel_time / all_pixel_times.size()) / thread_count;
  cout << "Sum of all pixel times: " << total_pixel_time << " s." << endl;
  cout << "Approximate Average Pixel Time (sum(all_pixel_times) / pixel_count / thread_count) = " << avg_pixel_time_approx << endl;
  cout << "Raytracing completed in " << raytrace_elapsed << " s." << endl;

  float sdev = 0;
  for(int i = 0 ; i < all_pixel_times.size() ; i++) {
	  sdev += (((all_pixel_times[i]/thread_count) - avg_pixel_time_approx) * 
      ((all_pixel_times[i]/thread_count) - avg_pixel_time_approx));
  }
  averageSamplingRate = double(total_samples)/double(total_pixels);
  float samples_sdev = 0;
  for(int i = 0 ; i < all_pixel_samples.size() ; i++) {
	  samples_sdev += ( (all_pixel_samples.at(i)-averageSamplingRate)*(all_pixel_samples.at(i)-averageSamplingRate) );
  }
  samples_sdev /= all_pixel_samples.size();
  samples_sdev = sqrt(samples_sdev);
  sdev /= all_pixel_times.size();
  sdev = sqrt(sdev);
  cout << "Pixel time sdev: " << sdev << endl;


  total_timer.endTimer();
  double total_time = total_timer.elapsedTime();
  file << "Raytrace time: " << raytrace_elapsed << endl;
  file << "Total_time: " << total_time << endl;
  file << "Time per pixel [end_time-start_time / pixel count]: " << time_per_pixel_overall << endl;
  file << "Time per pixel [(sum(pixel times) / pixel_count / thread_count): " << avg_pixel_time_approx << endl; 
  file << "Timer per pixel sdev: " << sdev << " s." << endl;
  file << "Scene polygon count: " << total_number_of_triangles << endl;
  file << "Polygons per ray: " << cis_stat_all.total_ray_poly_tests / total_rays << endl; //multiply by supersampling factor
  file << "Voxels per ray:" << cis_stat_all.total_voxels_traversed / total_rays << endl; //multiply by supersampling factor
  file << "Intersection tests: " << cis_stat_all.total_ray_poly_tests << endl;
  file << "Positive intersection test %: " << 100.0*double(cis_stat_all.found_poly_intersections)/double(cis_stat_all.total_ray_poly_tests) << "%." << endl;
	file << "Average Sampling Rate: " << averageSamplingRate << endl;
	file << "Sampling Rate Sdev: " << samples_sdev << endl;
  file << "Total samples: " << total_samples << "." << endl;
  file << "Average Ray Depth: " << double(total_depths) / double(total_samples) << endl;


  file.close();

  cout << "Process completed in " << total_time << " s." << endl;
  

  return MStatus::kSuccess;
}

//---------------------------------------------------------------------------
MString RayTracer::buildOutputDirectoryName()
{
  MString user_dir = MCommonSystemUtils::getEnv(MString("MAYA_APP_DIR"));
  MString img_dir = user_dir;
  MString last_char = img_dir.substring(img_dir.length()-1, img_dir.length()-1);
  if((MString("/") != last_char) && (MString("\\") != last_char)){
    img_dir += MString("/");
  }
  img_dir += MString("img/");
  return img_dir;
}

//---------------------------------------------------------------------------
bool RayTracer::intArg(const MArgList& args, unsigned int &indx, int & res)
{
  if (indx < args.length())
  {
    MStatus stat;
    indx++;
    res = args.asInt( indx, &stat );
    if (stat == MS::kSuccess){
      return true;
    }
  }
  return false;
}

//---------------------------------------------------------------------------
bool RayTracer::intArrayArg(const MArgList& args, unsigned int &indx, unsigned int count, std::vector<int> & res)
{
  int len = args.length();
  if (indx < len-1)
  {
    MStatus stat;
    for(int i=0; i<count; ++i){
      indx++;
      int res1 = args.asInt( indx, &stat );
      if (stat != MS::kSuccess){
        return false;
      }
      res.push_back(res1);
    }
    return true;
  }
  return false;
}


//---------------------------------------------------------------------------
bool RayTracer::fltArg(const MArgList& args, unsigned int &indx, float & res)
{
  if (indx < args.length())
  {
    MStatus stat;
    indx++;
    res = (float)args.asDouble( indx, &stat );
    if (stat == MS::kSuccess){
      return true;
    }
  }
  return false;
}

//---------------------------------------------------------------------------
bool RayTracer::strArg(const MArgList& args, unsigned int &indx, MString &res)
{
  if (indx < args.length())
  {
    MStatus stat;
    indx++;
    res = args.asString( indx, &stat );
    if (stat == MS::kSuccess){
      return true;
    }
  }
  return false;
}

//---------------------------------------------------------------------------
MStatus RayTracer::parseArgs(const MArgList& args,
                             int& width,
                             int& height,
                             int& voxelFactor,
                             int& superSample,
                             SuperSampleType& sst,
                             int& samplingRate,
                             float& tolerance,
                             int& misr,
                             int& masr,
                             int& rayDepth,
                             int& verbose, 
                             int& draw_x, 
                             int& draw_y, 
                             bool& drawOnly, 
                             int& threadCount)
{
	MStatus     	stat;
	MString     	arg;
	const MString	widthFlag				("-w");
	const MString	heightFlag			("-h");
	const MString	voxelFactorFlag	("-n");
	const MString	superSampleFlag	("-s");
  const MString superSampleTypeFlag ("-ss");
  const MString samplingRateFlag    ("-sr");
  const MString toleranceFlag       ("-t");
  const MString minSRFlag           ("-misr");
  const MString maxSRFlag           ("-masr");
  const MString rayDepthFlag        ("-rd");
	const MString	verboseFlag	        ("-v");
  const MString drawFlag            ("-draw");
  const MString drawOnlyFlag        ("-drawOnly");
  const MString threadCountFlag     ("-threads");


  width       = -1;
  height      = -1;
  voxelFactor = -1;
  superSample = -1;
  sst         = SuperSampleType_None;
  samplingRate = -1;
  tolerance   = -1.0f;
  misr        = -1;
  masr        = -1;
  rayDepth    = -1;
  verbose     = -1;
  draw_x      = -1;
  draw_y      = -1;
  drawOnly    = false;
  threadCount = -1;

	// Parse the arguments.
	for ( unsigned int i = 0; i < args.length(); i++ ) {
		arg = args.asString( i, &stat );
		if (!stat){
			continue;
    }

		if (arg == widthFlag)
			intArg(args, i, width);
		else if (arg == heightFlag)
			intArg(args, i, height);
		else if (arg == voxelFactorFlag)
			intArg(args, i, voxelFactor);
		else if (arg == superSampleFlag)
			intArg(args, i, superSample);
		else if (arg == superSampleTypeFlag) 
    {
      MString strSST;
			strArg(args, i, strSST);
      strSST.toLowerCase();
      if(SST_Uniform == strSST)
        sst = SuperSampleType_Uniform;
      else if(SST_Jittered == strSST)
        sst = SuperSampleType_Jittered;
      else if(SST_Random == strSST)
        sst = SuperSampleType_Random;
      else if(SST_Adaptive == strSST)
        sst = SuperSampleType_Adaptive;
      else
        sst = SuperSampleType_Unknown;
    }
		else if (arg == samplingRateFlag)
			intArg(args, i, samplingRate);
		else if (arg == toleranceFlag)
			fltArg(args, i, tolerance);
		else if (arg == minSRFlag)
			intArg(args, i, misr);
		else if (arg == maxSRFlag)
			intArg(args, i, masr);
		else if (arg == rayDepthFlag)
			intArg(args, i, rayDepth);
		else if (arg == verboseFlag)
			verbose = 1;
    else if (arg == drawFlag){
      std::vector<int> ints_out;
      intArrayArg(args, i, 2, ints_out);
      draw_x = ints_out[0];
      draw_y = ints_out[1];
      drawOnly = false;
    }
    else if (arg == drawOnlyFlag){
      std::vector<int> ints_out;
      intArrayArg(args, i, 2, ints_out);
      draw_x = ints_out[0];
      draw_y = ints_out[1];
      drawOnly = true;
    }
		else if (arg == threadCountFlag)
			intArg(args, i, threadCount);
		else {
			displayError(arg + ": unknown argument");
			//return MS::kFailure;
		}
	}
	return stat;
}

//---------------------------------------------------------------------------
bool RayTracer::getBoundingBox(const MDagPath& dagPath, bool quiet, BBox& bbox_)
{
  BBox bound_box = BBoxSentinel;

  MItMeshVertex it (dagPath);
  int count_vertices = 0;
  it.reset();
  while(!it.isDone())
  {
    MObject item = it.currentItem();
    MPoint pos = it.position(MSpace::kWorld);
    if(pos.x < bound_box.xmin) {
      bound_box.xmin = float(pos.x);
    }
    if(pos.x > bound_box.xmax) {
      bound_box.xmax = float(pos.x);
    }
    if(pos.y < bound_box.ymin) {
      bound_box.ymin = float(pos.y);
    }
    if(pos.y > bound_box.ymax) {
      bound_box.ymax = float(pos.y);
    }
    if(pos.z < bound_box.zmin) {
      bound_box.zmin = float(pos.z);
    }
    if(pos.z > bound_box.zmax) {
      bound_box.zmax = float(pos.z); 
    }
    it.next();
  }

  bbox_ = bound_box;
  return true;
}

//---------------------------------------------------------------------------
bool RayTracer::getBoundingBox(const std::vector<BBox>& bboxes, BBox& bbox_)
{
  if(bboxes.empty())
  { return false; }

  // loop through all bounding boxes
  BBox bound_box = bboxes[0];
  for(int i=1; i<bboxes.size(); ++i)
  {
    // compare each limit
    const BBox curr_bbox = bboxes[i];
    if(curr_bbox.xmin < bound_box.xmin) {
      bound_box.xmin = curr_bbox.xmin;
    }
    if(curr_bbox.xmax > bound_box.xmax) {
      bound_box.xmax = curr_bbox.xmax;
    }
    if(curr_bbox.ymin < bound_box.ymin) {
      bound_box.ymin = curr_bbox.ymin;
    }
    if(curr_bbox.ymax > bound_box.ymax) {
      bound_box.ymax = curr_bbox.ymax;
    }
    if(curr_bbox.zmin < bound_box.zmin) {
      bound_box.zmin = curr_bbox.zmin;
    }
    if(curr_bbox.zmax > bound_box.zmax) {
      bound_box.zmax = curr_bbox.zmax; 
    }
  }

  // send it out
  bbox_ = bound_box;
  return true;
}


//---------------------------------------------------------------------------
MColor RayTracer::calcColorAtPt(VoxelSpace& vxspc, float uv1, float uv2, 
                                const MFloatPoint& vertPt, const MFloatVector& vertNorm, const MFloatPoint& eyePoint, 
                                const Material& material, const std::vector<LightInfo>& light_info_list,
                                int mesh_idx, int face_idx, float bary1, float bary2,
                                CIStats& ci_stat,
                                const MColor& ref_color, const MColor& trans_color, float fres_kr)
{
  MColor ka_ = material.m_ambient;
  MColor kd_;
  if(material.m_isTextured){
    material.getBilinearFilteredPixelColor(uv1, uv2, kd_);
  } else {
    kd_ = material.m_diffuse;
  }

  //float diffuse_coeff = material.m_diffuseCoeff;
  //float reflectivity = material.m_reflectivity;
  float transparency = material.m_transparency.r;

  MColor ka = (1.0f-transparency) * ka_;          // (1-c)*ka_
  MColor kd = (1.0f-transparency) * kd_;          // (1-c)*kd_
  float fres_kt = transparency * (1.0f-fres_kr);  // c*kt_

  CIStats ci_stat_local;
  MColor pixel_color(0, 0, 0);
  for(int l=0; l<light_info_list.size(); ++l)
  {
    const LightInfo& li = light_info_list[l];
    MColor curr_component(0, 0, 0);
    switch (li.light_type)
    {
    case MFn::kAmbientLight:
      {
        // ka * ia
        curr_component = (ka * li.color) * li.intensity;
        break;
      }
    case MFn::kPointLight:
      {
        // test if this point visible from the light
        MFloatVector surf_light_vec = li.pos - vertPt;
        surf_light_vec.normalize();
        bool lightSeesPt = lightSeesPoint(vxspc, li.pos, -surf_light_vec, mesh_idx, face_idx, bary1, bary2, ci_stat_local);
        ci_stat += ci_stat_local;
				if(!lightSeesPt){
          continue;
        }

        // add the diffuse component kd*(N.L)*id
        float norm_vec_dot = surf_light_vec * vertNorm;
        norm_vec_dot = clampbot(norm_vec_dot, 0.0f);
        curr_component = (kd * norm_vec_dot * li.color);

        // add the specular component
        MFloatVector surf_light_ref(0.0f, 0.0f, 0.0f);
        if(norm_vec_dot>0.0f){
          surf_light_ref = 2*norm_vec_dot*vertNorm - surf_light_vec;
          surf_light_ref.normalize();

          MFloatVector surf_eye_vec = eyePoint - vertPt;
          surf_eye_vec.normalize();
          float norm_ref_dot = surf_light_ref*surf_eye_vec;
          norm_ref_dot = clampbot(norm_ref_dot, 0.0f);
          curr_component += (material.m_specular * pow(norm_ref_dot, material.m_specular.a) * li.color);
        }
        // multiply by the light intensity
        curr_component *= li.intensity;
        break;
      }
    case MFn::kDirectionalLight:
      {
				//MFloatVector vv(-1*li.dir.x,-1*li.dir.y,-li.dir.z);
        bool lightSeesPt = lightSeesPoint(vxspc, vertPt, li.dir, mesh_idx, face_idx, bary1, bary2, ci_stat_local, true);
        ci_stat += ci_stat_local;
        if(!lightSeesPt){
          continue;
        }

        // dir is the direction of the light. do the dot with -dir
        // the dot product becomes negative when incident ray and norm are in opposite dirs
        float norm_vec_dot = (-li.dir) * vertNorm;
        norm_vec_dot = clampbot(norm_vec_dot, 0.0f);
        curr_component = (kd * norm_vec_dot * li.color);

        // add the specular component
        MFloatVector surf_light_ref = 2*norm_vec_dot*vertNorm + li.dir;
        surf_light_ref.normalize();
        MFloatVector surf_eye_vec = eyePoint - vertPt;
        surf_eye_vec.normalize();
        float norm_ref_dot = surf_light_ref*surf_eye_vec;
        norm_ref_dot = clampbot(norm_ref_dot, 0.0f);
        curr_component += (material.m_specular * pow(norm_ref_dot, material.m_specular.a) * li.color);

        // multiply by the light intensity
        curr_component *= li.intensity;

        break;
      }
    }

    pixel_color += curr_component;
  }

  // add the reflective component
  MColor ref_comp = fres_kr*ref_color;
  pixel_color += ref_comp;

  // add the transmissive component
  MColor trans_comp = fres_kt*trans_color;
  pixel_color += trans_comp;

  // clamp to 1 to avoid overflows
  pixel_color.r = clamptop(pixel_color.r, 1.0f);
  pixel_color.g = clamptop(pixel_color.g, 1.0f);
  pixel_color.b = clamptop(pixel_color.b, 1.0f);


  return pixel_color;

}

//---------------------------------------------------------------------------
bool RayTracer::lightSeesPoint(const VoxelSpace& vxspc, const MFloatPoint& light_pos, const MVector& light_vec, 
                               int mesh_idx, int face_idx, float bary1, float bary2, 
                               CIStats& ci_stat, bool directional_light)
{
  ci_stat = CIStats();
	MFloatPoint hitPoint;
	//int hitFace;
  //int hitTriangle;
  int minMeshIdx=-1,minFaceIdx=-1;
	//float hitBary1;
  //float hitBary2;
  float hitBary3;
  //float hitRayParam;
  float minhitRayParam = FLT_MAX;
  float minBary1=-1.;
  float minBary2=-1.;

  bool intersects_any_mesh = false;
  MItDependencyNodes mesh_it(MFn::kMesh);
  int meshIdx = 0;

  bool bIntersectsOur = false;
  MFloatPoint pt;
  if(directional_light){
    pt.x = light_pos.x - vxspc.getDiagLength() * float(light_vec.x);
    pt.y = light_pos.y - vxspc.getDiagLength() * float(light_vec.y);
    pt.z = light_pos.z - vxspc.getDiagLength() * float(light_vec.z);
  }
  else {
    pt = light_pos;
  }
  triangle_object min_tri_obj;
  CIStats ci_stat_local;
  bIntersectsOur = vxspc.ourClosestIntersection(light_vec, pt, minhitRayParam, min_tri_obj, 
    hitPoint, minBary1, minBary2, hitBary3, 
    ci_stat_local, NULL);
  minMeshIdx = min_tri_obj.object;
  minFaceIdx = min_tri_obj.index;
  ci_stat += ci_stat_local;

  if(bIntersectsOur)
  {
    intersects_any_mesh = true;
    if(minMeshIdx != mesh_idx && directional_light)
      return false;
  }

  //if(directional_light && !intersects_any_mesh)
  //	return true;
  if(minMeshIdx != mesh_idx)
  { return false; }
  
  if(minFaceIdx != face_idx)
  { return false; }
  
  if(abs(minBary1-bary1)>0.0001 || abs(minBary2-bary2)>0.0001)
  { return false; }
  
  return true;
}

//---------------------------------------------------------------------------
void RayTracer::triangulateMesh(const MFnMesh& mesh)
{ 
    MString cmd("polyTriangulate -ch 0 "); 
    cmd += mesh.name(); 
    MGlobal::executeCommand( cmd ); 
    cmd = "select -cl";
    MGlobal::executeCommand( cmd );
}

//---------------------------------------------------------------------------
void RayTracer::triangulate_meshes()
{
  // loop through all the meshes
  MItDependencyNodes mesh_it(MFn::kMesh);
  for(; !mesh_it.isDone(); mesh_it.next())
  {
    // get the object
    MObject obj = mesh_it.item();
    MFnDagNode objNode(obj);
    cout << "triangulate_meshes(): " << objNode.name().asChar() << ": " <<  obj.apiTypeStr();

    // attach the function set
    MDagPath aPath = MDagPath::getAPathTo(obj);
    MFnMesh mesh(aPath);
    int pvc = mesh.polygonVertexCount(0);
    cout << " " << pvc << " vertices per polygon." << endl;

    // triangulate    
    triangulateMesh(mesh);
  }
}

//---------------------------------------------------------------------------
void RayTracer::get_scene_bounding_box(BBox& scene_bb, float scale_factor)
{
  MItDependencyNodes mesh_it(MFn::kMesh);
  std::vector<BBox> bbox_vector;
  BBox curr_bbox;
	for( ; !mesh_it.isDone(); mesh_it.next())
  {
    // get the object and a dagpath to it
    MObject obj = mesh_it.item();
    MFnDagNode objNode(obj);
    cout << "get_scene_bounding_box(): Including bounding box for: " << objNode.name().asChar() << " | " <<  obj.apiTypeStr() << endl;
    MDagPath aPath = MDagPath::getAPathTo(obj);

    // get and store the bounding box
    getBoundingBox(aPath, true, curr_bbox);
    bbox_vector.push_back(curr_bbox);
  }

  // calculate scene bounding box
  BBox all_obj_bbox;
  getBoundingBox(bbox_vector, all_obj_bbox);

  // expand it by a small factor
  float ctrx, ctry, ctrz;
  all_obj_bbox.center(ctrx,ctry,ctrz);
  float ext_x, ext_y, ext_z;
  all_obj_bbox.extents_half(ext_x, ext_y, ext_z);
  
  scene_bb.xmin = ctrx - (scale_factor*ext_x);
  scene_bb.xmax = ctrx + (scale_factor*ext_x);
  scene_bb.ymin = ctry - (scale_factor*ext_y);
  scene_bb.ymax = ctry + (scale_factor*ext_y);
  scene_bb.zmin = ctrz - (scale_factor*ext_z);
  scene_bb.zmax = ctrz + (scale_factor*ext_z);

}

//---------------------------------------------------------------------------
void RayTracer::create_meshinfo_list(std::vector<MeshInfo>& mesh_info_vec)
{
  MItDependencyNodes mesh_it(MFn::kMesh);
	for( ; !mesh_it.isDone(); mesh_it.next())
  {
    MeshInfo curr_mesh_geom;

    // get the object/name
    MObject obj = mesh_it.item();
    MFnDagNode objNode(obj);
    cout << "create_meshinfo_list(): creating MeshInfo struct for: " << objNode.name().asChar() << " | " <<  obj.apiTypeStr() << endl;
    curr_mesh_geom.m_name = objNode.name();
    
    // attach the function set
    MDagPath aPath = MDagPath::getAPathTo(obj);
    MFnMesh mesh(aPath);

    // store the dagpath
    curr_mesh_geom.m_dagPath = aPath;

    // m_pts on the whole mesh
    MFloatPointArray vertexList;
    mesh.getPoints(vertexList, MSpace::kWorld);
    for(unsigned int v=0; v<vertexList.length(); ++v){
      curr_mesh_geom.m_pts.push_back(vertexList[v]);
    }

    // per-polygon stuff
    MItMeshPolygon face_it(aPath);

    for(; !face_it.isDone(); face_it.next())
    {
      MeshInfo::FaceGeom curr_face_geom;

      // index
      int poly_idx = face_it.index();
      curr_face_geom.m_index = poly_idx;

      MIntArray polyVertList;
      mesh.getPolygonVertices(poly_idx, polyVertList);
      for(unsigned int pv=0; pv<polyVertList.length(); ++pv)
      {
        // point index in m_ptIdx
        int curr_pt_idx = polyVertList[pv];
        curr_face_geom.m_ptIdx.push_back(curr_pt_idx);

        // actual point (vertex) in m_pts
        MFloatPoint curr_pt = curr_mesh_geom.m_pts[curr_pt_idx];
        curr_face_geom.m_pts.push_back(curr_pt);

        // and m_normals
        MVector dcurr_normal;
        mesh.getFaceVertexNormal(poly_idx, polyVertList[pv], dcurr_normal, MSpace::kWorld);
        MFloatVector curr_normal(dcurr_normal);
        curr_face_geom.m_normals.push_back(curr_normal);
      }

      curr_mesh_geom.m_faceGeoms.push_back(curr_face_geom);
    }

    // push object on and then load the material 
    // (since we changed MImage to a pointer this might not be crucial)
    mesh_info_vec.push_back(curr_mesh_geom);
    mesh_info_vec[mesh_info_vec.size()-1].m_material.load(mesh);
  }

}

//---------------------------------------------------------------------------
void RayTracer::create_triangles_list(std::vector<std::vector<triangle>>& triangles_vector, int& triangle_count, bool verbose)
{
  MItDependencyNodes mesh_it(MFn::kMesh);
  triangles_vector.clear();

	size_t total_number_of_triangles = 0;  
	for( ; !mesh_it.isDone(); mesh_it.next())
  {
    // get the object
    MObject obj = mesh_it.item();
    MFnDagNode objNode(obj);
    cout << "create_triangles_list(): " << objNode.name().asChar() << " | " <<  obj.apiTypeStr() << endl;
    
    // attach the function set
    MDagPath aPath = MDagPath::getAPathTo(obj);
    MFnMesh mesh(aPath);

    // m_pts on the whole mesh
    MFloatPointArray vertexList;
    mesh.getPoints(vertexList, MSpace::kWorld);

    // save the vertices, triangles and normals for each mesh
    MIntArray Trianglesm,CountTriangles;
    mesh.getPoints(vertexList, MSpace::kWorld);
    int l = vertexList.length();
    mesh.getTriangles(CountTriangles, Trianglesm);
    int hn = Trianglesm.length();

    int ind = 0;
    std::vector<triangle> triangles;
    for(unsigned int t = 0 ; t < CountTriangles.length() ; t++)
    {
      int co = CountTriangles[t];
      for(int j = 0 ; j < co ; j++)
      {
        triangle tri;
        tri.p1 = vertexList[Trianglesm[ind*3]];
        tri.p2 = vertexList[Trianglesm[ind*3+1]];
        tri.p3 = vertexList[Trianglesm[ind*3+2]];
        triangles.push_back(tri);
        if(verbose){
          cout << "    poly " << t << " tri " << j << "/" << co << endl;
          cout << "         (" << tri.p1.x << ", " << tri.p1.y << ", " << tri.p1.z << ")" << endl;
          cout << "         (" << tri.p2.x << ", " << tri.p2.y << ", " << tri.p2.z << ")" << endl;
          cout << "         (" << tri.p3.x << ", " << tri.p3.y << ", " << tri.p3.z << ")" << endl;
        }
        ind++;
      }
    }
    total_number_of_triangles += triangles.size();
    triangles_vector.push_back(triangles);
  }

  triangle_count = (int)total_number_of_triangles;
}

//---------------------------------------------------------------------------
MStatus RayTracer::get_renderplane_info(int width, int height, RenderPlaneInfo& rp_info)
{

  // look at all cameras
  bool found_non_default_cam = false;
  int ImgResX = width;
  int ImgResY = height;
  MPoint LeftBottom;
  double PixelHeight;
  double PixelWidth;
  MFloatPoint eyePoint;
  MVector viewDirection;
  MVector upDirection;
  MVector rightDirection;
  bool isOrtho = false;

  double ImgPlaneH;
  double ImgPlaneW;
  double aspectRatio;

  double left;
  double right;
  double bottom;
  double top;
  double windowAspect = 1.0;
  
  MPoint imagePlaneCenter;

  MItDependencyNodes cam_it(MFn::kCamera);
  while(!cam_it.isDone())
  {
    MObject obj = cam_it.item();
    MDagPath aPath = MDagPath::getAPathTo(obj);
    MStatus status;
    MFnCamera camera (aPath, &status);
    if ( !status ) {
      status.perror("MFnCamera constructor");
      continue;
    }
    std::cout << "Found camera " << camera.name().asChar();

    // skip if it's one of the default cameras
    if((MString("perspShape") == camera.name()) ||
      (MString("topShape") == camera.name()) ||
      (MString("frontShape") == camera.name()) ||
      (MString("sideShape") == camera.name()))
    {
      std::cout << "[default camera; skipping...]" << std::endl;
      cam_it.next();
      continue;
    }

    std::cout << std::endl;
    found_non_default_cam = true;
    camera.getRenderingFrustum(windowAspect, left, right, bottom, top);
    MPoint dEyePoint = camera.eyePoint(MSpace::kWorld);
    eyePoint.setCast(dEyePoint);

    viewDirection =	camera.viewDirection(MSpace::kWorld); 
    viewDirection.normalize();
    upDirection = camera.upDirection(MSpace::kWorld);
    upDirection.normalize();
    rightDirection = camera.rightDirection(MSpace::kWorld);
    rightDirection.normalize();

    double focalLengthCM = camera.focalLength() / 10.0;
    // film aperture is in inches.
    double horizontalFilmApertureCM = camera.horizontalFilmAperture() * 2.54; 
    aspectRatio = camera.aspectRatio();

    ImgPlaneW = horizontalFilmApertureCM;
    ImgPlaneH = ImgPlaneW / aspectRatio;

    if(camera.isOrtho())
    {
      isOrtho = true;
      float orthographic_width = (float)camera.orthoWidth();
      aspectRatio = camera.aspectRatio();
      ImgPlaneW = orthographic_width;
      ImgPlaneH = ImgPlaneW / aspectRatio;
    }

    //viewAngle = 2*atan(horizontalFilmAperture/(2*focalLength));
    //H = 2*focalLength*tan(viewAngle/2);

    //W = H*aspectRatio;
    MVector tmp = viewDirection*focalLengthCM;
    MPoint p(eyePoint);
    p.x = p.x + tmp.x; p.y = p.y + tmp.y; p.z = p.z + tmp.z;
    imagePlaneCenter = p;

    LeftBottom = imagePlaneCenter - (ImgPlaneW/2)*rightDirection;
    LeftBottom = LeftBottom - (ImgPlaneH/2)*upDirection;
    PixelWidth = ImgPlaneW/ImgResX;
    PixelHeight = ImgPlaneH/ImgResY;

    break;
  }

  if(!found_non_default_cam)
  {
    std::cout << "Found no non-default cameras. Exiting..." << std::endl;
    return MStatus::kFailure;
  }

  rp_info.ImgResX = width;
  rp_info.ImgResY = height;
  rp_info.LeftBottom = LeftBottom;
  rp_info.PixelHeight = PixelHeight;;
  rp_info.PixelWidth = PixelWidth;
  rp_info.eyePoint = eyePoint;
  rp_info.viewDirection = viewDirection;
  rp_info.upDirection = upDirection;
  rp_info.rightDirection = rightDirection;
  rp_info.isOrtho = isOrtho;
  rp_info.ImgPlaneH = ImgPlaneH;
  rp_info.ImgPlaneW = ImgPlaneW;
  rp_info.aspectRatio = aspectRatio;


  std::cout << (isOrtho?"Orthographic":"Perspective") << " Cam Properties: " << std::endl;
  std::cout << "  {Left, Right, Bot, Top} = {" << left << ", " << right << ", " << 
    bottom << ", " << top << "}" << std::endl;


  return MStatus::kSuccess;
}

//---------------------------------------------------------------------------
void RayTracer::create_lightinfo_list(std::vector<LightInfo>& light_info_list)
{
  light_info_list.clear();
  MItDependencyNodes light_it(MFn::kLight);
  while(!light_it.isDone())
  {
    MDagPath light_path = MDagPath::getAPathTo(light_it.item());
    MStatus stat;
    MFnLight curr_light(light_path, &stat);
    if(!stat)  {
      std::cout << "Not a light." << std::endl;
      continue;
    }

    // is a light, continue.
    std::cout << "create_lightinfo_list(): Found light of type " << light_it.item().apiTypeStr() << " ";

    MFn::Type light_api_type = light_it.item().apiType();
    if((MFn::kAmbientLight != light_api_type) && 
      (MFn::kDirectionalLight != light_api_type) && 
      (MFn::kPointLight != light_api_type))
    { 
      std::cout << "[Not Ambient, Directional or Point, skipping...]" << std::endl;
      continue;
    }

    MColor light_color = curr_light.color();
    std::cout << "(R,G,B) = (" << light_color.r << ", " << light_color.g << ", "
      << light_color.b << ")" << std::endl;

    MMatrix light_mat = light_path.inclusiveMatrix(&stat);
    if(!stat)
    {
      std::cout << "Error getting inclusive matrix." << std::endl;
    }

    MPoint origin(0, 0, 0, 1);
    MVector xaxis(1, 0, 0);
    MVector yaxis(0, 1, 0);
    MVector zaxis(0, 0, 1);
    MVector origin_world = origin * light_mat;
    MVector xaxis_world = xaxis * light_mat;
    MVector yaxis_world = yaxis * light_mat;
    MVector zaxis_world = zaxis * light_mat; 

    // to find the direction the light is pointing.
    MVector negZaxis(0, 0, -1);
    MVector negZaxis_world = negZaxis * light_mat;

    MFnDagNode transform(light_path.transform());
    MMatrix transformMat = transform.transformationMatrix();
    MTransformationMatrix light_transform(transformMat);
    MFloatPoint  light_pos = light_transform.getTranslation(MSpace::kWorld);

    LightInfo curr_light_info;
    curr_light_info.light_path = light_path;
    curr_light_info.light_type = light_api_type;
    curr_light_info.pos = light_pos;
    curr_light_info.dir = MFloatVector(0, 0, 0);
    curr_light_info.color = curr_light.color();
    curr_light_info.intensity = curr_light.intensity();
    if(MFn::kDirectionalLight == light_api_type){
      curr_light_info.dir = negZaxis_world;
    }

    light_info_list.push_back(curr_light_info);
    light_it.next();
  }
}

void RayTracer::computeChiInvForBetas()
{
	m_chiInv_R = new double* [CHI_INV_COUNT];
	int c = 0;
	for(double beta = 1.0/double(CHI_INV_COUNT) ; beta < 1 ; beta += 1.0/double(CHI_INV_COUNT))
	{
		m_chiInv_R[c] = new double [m_arg_maxSR+1];
		for(int i = 1 ; i <= m_arg_maxSR ; i++)
		{
			m_chiInv_R[c][i] = inverseChi(i,beta);
		}
		c++;
	}
}

//---------------------------------------------------------------------------
//void RayTracer::cleanUpChiInv()
//{
//  if(NULL != m_chiInv_R){
//    for(int i=0; i<CHI_INV_COUNT; ++i){
//      delete [] m_chiInv_R[i];
//      m_chiInv_R[i] = NULL;
//    }
//    delete [] m_chiInv_R;
//    m_chiInv_R = NULL;
//  }
//  if(NULL != m_chiInv_G){
//    for(int i=0; i<CHI_INV_COUNT; ++i){
//      delete [] m_chiInv_G[i];
//      m_chiInv_G[i] = NULL;
//    }
//    delete [] m_chiInv_G;
//    m_chiInv_G = NULL;
//  }
//  if(NULL != m_chiInv_B){
//    for(int i=0; i<CHI_INV_COUNT; ++i){
//      delete [] m_chiInv_B[i];
//      m_chiInv_B[i] = NULL;
//    }
//    delete [] m_chiInv_B;
//    m_chiInv_B = NULL;
//  }
//}

//---------------------------------------------------------------------------
bool RayTracer::RayTrace(const MFloatPoint& raySource, const MFloatVector& rayDir, int traceDepth, float curr_refraction_idx,
                         MColor& out_color, 
                         CIStats& ci_stats_, int& total_rays_, std::vector<int>& sample_ray_depths,
                         bool draw, triangle_object* tri_skip)
{
  ci_stats_ = CIStats();
  total_rays_ = 0;

  // set up spaces
  char spc[30];
  if(draw){
    memset(spc, (int)' ', traceDepth+1);
    spc[traceDepth+1] = '\0';
  }

  if(traceDepth > m_arg_rayDepth)
  {
    if(draw){
      cout << spc << "RayTrace(td=" << traceDepth << "): Max depth exceeded." << endl;
    }
    out_color = MColor(0.0f, 0.0f, 0.0f);
    int depth_to_return = traceDepth-1;
    sample_ray_depths.push_back(depth_to_return);
    return true;
  }



  if(draw){
    // box at origin
    RayTracer::s_viz.drawCube(raySource.x, raySource.y, raySource.z, 0.1, 0.1, 0.1);
    RayTracer::s_viz.deselectAll();
    cout << spc << "RayTrace(td=" << traceDepth << "): SRC: " << raySource.x << "," << raySource.y << "," << raySource.z << endl;
  }


  MFnMesh minMesh;
  bool bIntersectedAtAll = false;
  bool bIntersectedAtAllCheck = false;

  float minHitRayParamOur, BaryOur1, BaryOur2, BaryOur3;
  int minMeshOur, minHitTriangleOur; 
  MFloatPoint minHitPointOur;
  float2 uv1,uv2,uv3,minuv;

  triangle_object min_tri;
  CIStats ci_stats_local;
  
  bool bRayHit = m_vxspc.ourClosestIntersection(rayDir, raySource, minHitRayParamOur, min_tri, minHitPointOur, 
    BaryOur1, BaryOur2, BaryOur3, ci_stats_local, tri_skip);
  minMeshOur = min_tri.object;
  minHitTriangleOur = min_tri.index;
  ci_stats_ += ci_stats_local;
  ++total_rays_;

  if(!bRayHit)
  {
    if(draw){
      MPoint origin(raySource);
      MPoint dest = origin + 0.5*MVector(rayDir.x,rayDir.y,rayDir.z);
      RayTracer::s_viz.drawRay(origin, dest);
      RayTracer::s_viz.deselectAll();
      cout << spc << "RayTrace(td=" << traceDepth << "): (NO HIT)" << endl;
    }

    // return black
    out_color = MColor(0.0f, 0.0f, 0.0f);
    sample_ray_depths.push_back(traceDepth);
    return true;
  }

  // if we're here, we intersected something

  // if we're debugging closestIntersection, return white at the top level
  static int debugging_closest_intersection = 0;
  if(debugging_closest_intersection && (0 == traceDepth))
  {
    // return white
    out_color = MColor(1.0f, 1.0f, 1.0f);
    sample_ray_depths.push_back(traceDepth);
    return true;
  }

  MPoint vtx1 = m_mesh_info_vec[minMeshOur].m_faceGeoms[minHitTriangleOur].m_pts[0];
  MPoint vtx2 = m_mesh_info_vec[minMeshOur].m_faceGeoms[minHitTriangleOur].m_pts[1];
  MPoint vtx3 = m_mesh_info_vec[minMeshOur].m_faceGeoms[minHitTriangleOur].m_pts[2];
  const MString str = (m_mesh_info_vec[minMeshOur].m_material.set);

  // overall mesh info struct & material & face geom
  const MeshInfo& curr_mesh_info = m_mesh_info_vec[minMeshOur];
  const Material& curr_mat = curr_mesh_info.m_material;
  const MeshInfo::FaceGeom& curr_face_geom = curr_mesh_info.m_faceGeoms[minHitTriangleOur];
  const MDagPath& curr_mesh_dagpath = curr_mesh_info.m_dagPath;

  MFnMesh curr_mesh(curr_mesh_dagpath);
  MStatus s1 = curr_mesh.getPolygonUV(minHitTriangleOur, 0, uv1[0], uv1[1], &str);
  MStatus s2 = curr_mesh.getPolygonUV(minHitTriangleOur, 1, uv2[0], uv2[1], &str);
  MStatus s3 = curr_mesh.getPolygonUV(minHitTriangleOur, 2, uv3[0], uv3[1], &str);
   
  minuv[0] = BaryOur1*uv1[0] + BaryOur2*uv2[0] + (1.0f - BaryOur1 - BaryOur2)*uv3[0];
  minuv[1] = BaryOur1*uv1[1] + BaryOur2*uv2[1] + (1.0f - BaryOur1 - BaryOur2)*uv3[1];

  MFloatPoint triV1 = curr_face_geom.m_pts[0];
  MFloatPoint triV2 = curr_face_geom.m_pts[1];
  MFloatPoint triV3 = curr_face_geom.m_pts[2];
  MFloatVector normV1 = curr_face_geom.m_normals[0];
  MFloatVector normV2 = curr_face_geom.m_normals[1];
  MFloatVector normV3 = curr_face_geom.m_normals[2];

  MFloatPoint ptWeighted = (triV1*BaryOur1 + triV2*BaryOur2 + triV3*(1.0f - BaryOur1 - BaryOur2));
  MFloatVector normWeighted = (BaryOur1*normV1 + BaryOur2*normV2 + (1.0f - BaryOur1 - BaryOur2)*normV3);
  normWeighted.normalize();

  if(draw){
    MPoint origin(raySource);
    MPoint dest(ptWeighted);
    RayTracer::s_viz.drawRay(origin, dest);
    RayTracer::s_viz.deselectAll();
    RayTracer::s_viz.drawCube(dest.x, dest.y, dest.z, 0.1, 0.1, 0.1);
    RayTracer::s_viz.deselectAll();
    cout << spc << "RayTrace(td=" << traceDepth << "): DST: " << dest.x << "," << dest.y << "," << dest.z << endl;
  }

  // figure out whether we're entering or exiting the material
  bool ray_entering_material = true;
  float incident_dot_norm = rayDir * normWeighted;
  if(incident_dot_norm > 0.0f)
  {
    // if incident and normal point in same dir, incident is leaving the mesh
    ray_entering_material = false;
  }

  // eta_i is the refraction index of the current medium
  float eta_i = curr_refraction_idx;

  // if rayDir is same dir as faceNorm, transmitted ray is exiting the object, so eta_t is 1
  // if transmitted ray is entering the object, eta_t takes its refraction index (otherwise, 
  // it's exiting the material, in which case it's 1.0f)
  float eta_t;
  if(ray_entering_material){
    eta_t = curr_mesh_info.m_material.m_refractiveIndex;
  } else {
    eta_t = 1.0f;
  }

  if(draw){
    cout << spc << "RayTrace(td=" << traceDepth << "): eta_i: " << eta_i << endl;
    cout << spc << "RayTrace(td=" << traceDepth << "): eta_t: " << eta_t << endl;
  }

  MFloatVector normToUse = normWeighted;
  if(!ray_entering_material){
    // if ray is exiting the material, flip the normal
    normToUse = -normWeighted;
  }


  // calculate reflection -> refl_dir
  MFloatVector refl_dir;
  refl_dir = reflection(-rayDir, normToUse);  

  if(draw){
    cout << spc << "RayTrace(td=" << traceDepth << "): REFL: " << refl_dir.x << "," << refl_dir.y << "," << refl_dir.z << endl;
  }
  
  // calculate transmission-> trans_dir
  float theta_i;
  MFloatVector trans_dir;
  bool valid_result;
  trans_dir = transmission(-rayDir, normToUse, eta_i, eta_t, valid_result, &theta_i);
  if(draw){
    if(valid_result){
      cout << spc << "RayTrace(td=" << traceDepth << "): TRNS1: " << trans_dir.x << "," << trans_dir.y << "," << trans_dir.z <<
        " theta_i: " << toDegrees(theta_i) << endl;
    } else {
      cout << spc << "RayTrace(td=" << traceDepth << "): TRNS1: (critical angle)" << " theta_i: " << toDegrees(theta_i) << endl;
    }
  }

  // calculate fresnel parameters
  float fres_kr0 = fresnel_kr0(eta_i, eta_t);
  float fres_kr = fresnel_kr(fres_kr0, theta_i);
  if(!valid_result){
    fres_kr = 1.0f;
  }
  if(draw){
    cout << spc << "RayTrace(td=" << traceDepth << "): kr0: " << fres_kr0 << endl;
    cout << spc << "RayTrace(td=" << traceDepth << "): kr(" << toDegrees(theta_i) << "): " << fres_kr << endl;
    if(!valid_result){
      cout << spc << " (Critical angle->total internal reflection)" << endl;
    }
  }


  // raytrace for reflected color. eta_i stays the same since it's in the same medium
  MColor refl_color;
  int ray_count_local;
  if(draw){
    cout << spc << "RayTrace(td=" << traceDepth << "): Starting Reflective RayTrace" << endl;
  }
  RayTrace(ptWeighted, refl_dir, traceDepth+1, eta_i, refl_color,
    ci_stats_local, ray_count_local, sample_ray_depths,
    draw, &min_tri);

  ci_stats_ += ci_stats_local;
  total_rays_ += ray_count_local;

  // raytrace for transmitted color. eta_t is the ref_idx of the new medium
  MColor trans_color;
  if(NearlyEq(curr_mat.m_transparency.r, 0.0f) && 
    NearlyEq(curr_mat.m_transparency.g, 0.0f) &&
    NearlyEq(curr_mat.m_transparency.r, 0.0f))
  {
    // if it's not transparent at all, no reason to do a RayTrace for transmission
    trans_color = MColor(0,0,0);
  } 
  else 
  {
    if(valid_result){
      if(draw){
        cout << spc << "RayTrace(td=" << traceDepth << "): Starting Transmissive RayTrace" << endl;
      }
      RayTrace(ptWeighted, trans_dir, traceDepth+1, eta_t, trans_color,
        ci_stats_local, ray_count_local, sample_ray_depths, draw, &min_tri);
      ci_stats_ += ci_stats_local;
      total_rays_ += ray_count_local;
    } else {
      if(draw){
        cout << spc << "RayTrace(td=" << traceDepth << "): transmission invalid. Color = (0,0,0)" << endl;
      }
      trans_color = MColor(0,0,0);
    }
  }  

  // Combine colors
  // phong (per-pixel) of local surface
  MColor pixel_color = RayTracer::calcColorAtPt(m_vxspc, minuv[0], minuv[1], ptWeighted, normToUse, raySource, curr_mat, m_light_info_list,
    minMeshOur, minHitTriangleOur, BaryOur1, BaryOur2, ci_stats_local, refl_color, trans_color, fres_kr);
  ci_stats_ += ci_stats_local;

  // return the pixel color
  out_color = pixel_color;
  return true;
}